!
!  QUICKSORT - recursive version of the QuickSort algorithm 
!
!  by  Wirth,N: Algorithm + Data Structure = Programs, Prentice-Hall, 1975
!
!  Copyright © 1997-2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.


module mquick

  public :: quicksort, quicksorti
  private :: trid, tridi

contains

  subroutine quicksort(a)

    real, dimension(:), intent(in out) :: a

    call trid(1,size(a),a)

  end subroutine quicksort



  recursive subroutine trid(l,r,a)

    integer,intent(in) :: l,r
    real, dimension(:),intent(in out) :: a

    integer i,j
    real x,w     ! internal buffers

    i = l
    j = r
    x = a((l+r)/2)

    do
       do while ( a(i) < x )
          i = i + 1
       end do

       do while ( x < a(j) )
          j = j - 1
       end do

       if( i <= j )then
          w = a(i); a(i) = a(j); a(j) = w
          i = i + 1; j = j - 1
       endif

       if( i > j ) exit
    end do

    if( l < j ) call trid(l,j,a)
    if( i < r ) call trid(i,r,a)

  end subroutine trid

!--------------------------------------------------------------------------

  subroutine quicksorti(a,idx)
    
    real, dimension(:), intent (in out) :: a
    integer, dimension(:), intent (in out) :: idx
    integer :: n

    n = size(a)

    if( n /= size(idx) ) stop 99

    if( n > 1 ) call tridi(1,n,a,idx)

  end subroutine quicksorti


  recursive subroutine tridi(l,r,a,idx)

    integer,intent( in ) :: l, r 
    integer, dimension(:), intent(in out) :: idx
    real, dimension(:), intent(in out) :: a

    integer i,j,ii
    real x,w     ! internal buffers

    i = l; j = r
    x = a((l+r)/2)

    do
       do while ( a(i) < x )
          i = i + 1
       end do

       do while ( x < a(j) )
          j = j - 1
       end do

       if( i <= j )then
          w = a(i); a(i) = a(j); a(j) = w
          ii = idx(i); idx(i) = idx(j); idx(j) = ii
          i = i + 1; j = j - 1
       endif

       if( i > j ) exit
    end do

    if( l < j ) call tridi(l,j,a,idx)
    if( i < r ) call tridi(i,r,a,idx)

  end subroutine tridi


end module mquick

